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FOREWORD 


This  report  is  the  result  of  research  performed  in  the 
Computational  Aerodynamics  Group,  Aerodynamics  and  Airframe  Branch, 
Aeromechanics  Division,  Flight  Dynamics  Laboratory  from  August  1983  to 
Noventer  1984.  This  report  was  prepared  by  Dr.  Miguel  R.  Visbal, 
Visiting  Scientist,  under  work  unit  2307  N603,  "Computational  Fluid 
Dynamics,"  with  Dr.  Wylbur  Hankey  as  the  Task  Manager.  The  report  was 
submitted  in  March  1985. 

The  author  would  like  to  express  his  appreciation  to 
Dr.  Wylbur  Hankey  and  Dr.  Joseph  Shang  for  their  valuable  technical 
guidance  and  advice  during  this  work.  My  association  with  all  the 
meiibers  of  the  Computational  Aerodynamics  Group  has  been  both  personally 
and  technically  rewarding. 

Computer  time  for  the  present  study  was  provided  in  part  by  the 
NASA  Ames  Research  Center. 
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SECTION  I 


INTRODUCTION 

In  recent  years,  the  numerical  simulation  of  high  Reynolds  number 
transonic  flows  over  airfoils  has  been  presented  in  numerous 

publications.  A  comprehensive  review  of  the  subject  can  be  found,  for 
example,  in  References  1  and  2,  and  therefore  is  not  repeated  here. 

This  earlier  work  clearly  illustrates  the  potential  of  Navier-Stokes 
numerical  methods  for  evolving  into  a  reliable  predictive  design 

technique.  However,  for  this  goal  to  be  attained,  two  major  problem 
areas  still  require  further  investigation.  First,  an  assessment  of  the 
accuracy  of  the  computed  solutions  must  be  conducted  in  order  to  clarify 
the  uncertainties  due  to  grid  resolution,  boundary  condition  formulation 
and  numerical  damping.  Second,  a  suitable  turbulence  model  must  be 

developed,  capable  of  accurately  describing  flows  containing  transition. 
Separation,  wakes  and  shock  wave/boundary  layer  interactions.  This 
second  problem  area  clearly  constitutes  the  pacing  item  in  computational 
aerodynamics  and  will  most  likely  require  a  continued  process  of 
numerical  evaluation  of  different  turbulence  modeling  formulations.  For 
this  evaluation  process  to  be  meaningful,  the  degree  of  uncertainty 
associated  with  the  numerical  simulations  must  be  first  established.  It 
should  be  noted,  however,  that  numerical  resolution  and  turbulence 
modeling  cannot  be  entirely  separated  since  grid  resolution  is  dependent 
on  the  turbulence  model  employed. 


Since  "exact"  solutions  are  not  generally  available  for  complex 
transonic  flows  of  interest,  grid  refinement  studies  constitute  the  only 
suitable  alternative  for  accuracy  assessment.  Even  for  two-dimensional 
flows,  this  approach  can  be  very  costly  in  terms  of  computer  resources 
and  therefore  the  accuracy  of  the  numerical  solutions  is  customarily 
judged  by  comparison  with  available  experimental  data.  However, 
comparing  computed  and  experimental  results  can  be  inconclusive  due  to 
(1)  lack  of  numerical  resolution,  (2)  uncertainties  in  the  measurements 
(e.g.  wind  tunnel  and  probe  interference),  and  (3)  the  inability  to 
isolate  numerical  errors  from  those  due  to  turbulence  modeling.  With  the 
advent  of  powerful  computers,  more  detailed  calculations  of 

3 

two-dimensional  transonic  viscous  flows  are  now  possible  ,  and  should 
help  clarify  some  of  the  uncertainties  inherent  to  the  numerical 
simulations. 

With  the  above  background  as  motivation,  the  present  research 
critically  examines  several  aspects  of  the  numerical  solution  of  the 
Navier-Stokes  equations  for  high  Reynolds  number  transonic  airfoil 

flows.  This  problem  embodies  many  interesting  flow  features  such  as 

leading  and  trailing  edge  regions,  wake  and  shock/boundary  layer 

interactions.  The  particular  configuration  considered  was  a  modified 

4 

Whitcomb  supercritical  airfoil,  denoted  as  DSMA  523.  This  airfoil 
exhibits  significant  rear-loading  and  a  strong  viscous-inviscid 
interaction  as  compared  with  more  conventional  sections  having  little  or 
no  aft-camber.  In  addition,  this  airfoil  was  selected  as  one  of  the  test 

5 

cases  for  the  1980  Stanford  conference  on  complex  turbulent  flows. 


Computations  were  performed  using  the  mass-averaged  Navier-Stokes 

equations^  expressed  in  terms  of  general  curvilinear  coordinates  and 

with  turbulence  represented  by  an  algebraic  eddy  viscosity  model. ^  The 

0 

governing  equations  were  solved  in  nearly -orthogonal  body-fitted  grids 

9 

employing  MacCormack's  explicit  scheme  and  the  implicit  Beam-Warming 
algorithm. These  two  most  commonly  used  numerical  schemes  were 
chosen  for  their  intrinsic  differences  in  solving  the  system  of  governing 
equations.  A  direct  comparison  between  the  two  algorithms  was  performed 
under  identical  mesh  systems,  boundary  conditions  and  turbulence  model. 
Such  a  one-to-one  comparison,  ^e  author  believes,  has  not  been 
previously  documented  in  the  literature. 

The  main  objectives  in  this  investigation  of  transonic  airfoil 
flows  can  be  summarized  as  follows: 

(1)  A  comparative  study  of  the  relative  accuracy  between  an 
explicit  and  an  implicit  Mavier-Stokes  code. 

(2)  Effects  of  grid  refinement,  numerical  damping  and  farfield 
boundary  placement  on  the  computed  flowfields. 

(3)  Comparison  of  the  results  obtained  with  the  Euler,  the 
thin-layer  and  the  full  Navier-Stokes  equations. 

(4)  Evaluation  of  the  algebraic  turbulent  eddy  viscosity  model  by 
detailed  comparison  of  computed  and  experimental  data  for  both 
subcritical  and  supercritical  flows. 


SECTION  II 


MATHEMATICAL  DESCRIPTION  OF  THE  FLOW 

In  this  investigation,  viscous  transonic  airfoil  flows  are 

simulated  by  means  of  the  two-dimensional  compressible  Navier-Stokes 
equations.  The  governing  equations  are  written  in  terms  of  general 
curvilinear  coordinates  and  mass-averaged  variables,  with  turbulence 
incorporated  through  an  algebraic  eddy  viscosity  model.  The  particular 
form  of  the  flow  equations,  turbulence  model  and  boundary  conditions 

employed  are  presented  in  this  section. 

1 .  GOVERNING  EQUATIONS 

The  governing  equations  are  taken  to  be  the  two-dimensional 

compressible  Navier-Stokes  equations  expressed  in  terms  of  mass-averaged 

variables®  and  general  curvilinear  coordinates  with  turbulence 

represented  by  an  algebraic  eddy  viscosity  model.  Two  different  forms  of 

the  equations  are  employed  in  this  research.  The  explicit  numerical 

algorithm  (described  in  Section  III)  utilizes  the  chain-rule  conservative 

form^\  while  the  implicit  scheme  solves  the  strong  conservative 

12 

formulati  on. 

In  terms  of  general  curvilinear  coordinates  (E.n),  the  chain-rule 
conservative  Navier-Stokes  equations  can  be  written  as  follows: 


where 


(2.2) 


pU2+p-T^^ 

pUV-Txy 

(pe  +  p)u 


pv^  +  P  - 
(pe  +  p)v  - 


(2.3) 


W  2(P  +  e)u^  +  A_^(u^+  Vy) 

T^y  =  (P  +  e)(Uy  +  v^) 

tww  =  2(p  +  e)v^,  +  A_(u^  +  v) 


(2.4) 


UT 

+ 

XX 

xy 

UT 

+ 

VT„., 

xy 

yy 

rt 


yy  p  ' 


The  variables  u  and  v  denote  respectively  the  x  and  y  velocity 
components.  The  density  p  ,  static  pressure  p  and  the  absolute 
temperature  T  satisfy  the  equation  of  state  for  a  perfect  gas 


p  =  pRT  =  (y  -  1)  [pe  -  p  (u^  +  v^)]  (2.5) 

where  e  is  the  total  energy  per  unit  mass,  R  is  the  gas  constant  (1716 
ft /sec  -  "R  for  air)  and  Y  (=1.4)  is  the  ratio  of  specific  heats. 
The  dynamic  molecular  viscosity  y  is  obtained  from  Sutherland's  formula. 
The  turbulent  eddy  viscosity  e  is  given  by  the  eddy  viscosity  model  which 
will  be  discussed  later.  Stokes  hypothesis  is  assumed  (i.e.,X.|.  =  -2/3 
(y+e  )).  The  molecular  Prandtl  number  Pr,  the  turbulent  Prandtl  number 


o 


Pr^  and  the  specific  heat  Cp  are  taken  to  be  constant  (Pr  =  0.72, 

Pr.  =  0.9  and  c„  =  6006  ft^/sec^  -  *R) 
t  p 

The  quantities  C  C  n  n  denote  the  transformation 

X ,  y  X  y 

metrics  (see  Figure  1  for  an  sketch  of  the  basic  transformation) 


C  =  y  J 


"x  ■ 


"y  '  ’'s'' 


(2.6) 


where  J  is  the  transformation  Jacobian 


J  = 


Mi 


¥(x 


(2.7) 


The  Navier-Stokes  equations  (2.1)  can  be  cast  in  strong 
12 

conservation  form  as  follows 


93^  +  ^  ^  n 

3t  3n  ■ 

where  a 

f  =  (^x^"  ^  ^y^)/J  (2.9) 

G  =  (n^F  +  qyG)/J 


In  order  to  facilitate  the  implementation  of  the  implicit  algorithm 


(Section  III),  Equation  (2.8)  is  rewritten  in  the  following  fashion 


11,1'  denote  the  contravar iant  velocities 


m 


U  -  +  CyV 

V  =  n^u  +  n, V 

X  y 


(2.16) 


and  the  viscous  coefficients  b^. ,  c.,  d.  are 


l>,  =  (u  ♦  C)(fcj  +  5>) 

K,  '  j<U  + 

l>,  -  (U  ♦  e)(cj  ♦ 

\  *  f: 


=  -(y  +  e)(|-  +  Cyiiy) 

C^  =  -(y  +  £)(■!  -  CyH^) 

S  =  (U  *  e)(My  -  |Vx> 

c,  =  -(l'  *^)(5A 

S  ^  Vy) 

^  '"t 


=  (P  +  e)(3  +  n^) 

=  +  e)n^ny 

=  (y  +  G)(n^  +  y  np 


(2.17) 


(2.18) 


(2.19) 
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By  neglecting  the  viscous  terms  in  the  direction  parallel  to  a  solid 
surface,  the  so  called  thin-layer  approximation^  is  obtained.  Assuming 
a  body -conforming  curvilinear  coordinate  system,  with  the  C  and  n 
directions  defined  as  shown  in  Figure  1,  the  thin-layer  Navier-Stokes 
equations  become 


3q  ^  9Ei  3^2 
9t  3^  3n 


3W  2  /  ^  \ 


(2.20) 


Finally,  the  Euler  equations  are  obtained  by  neglecting  all  viscous 
terms  (i.e.,  by  setting  the  right-hand-side  of  Equation  (2.20)  equal  to 
zero) . 

2.  BOUNDARY  AND  INITIAL  CONDITIONS 


In  order  to  completely  define  the  problem,  suitable  boundary  and 
initial  conditions  must  be  specified.  Referring  to  the  airfoil 
computational  domain  shown  in  Figure  1,  the  following  boundary  conditions 
are  prescribed. 

Along  the  outer  boundary  ABC,  freestream  conditions  are  given  for 
all  the  flow  variables.  On  the  downstream  boundaries  AD  and  HC, 


fi] 


Since  the  above  farfield  conditions  are  only  approximate,  the  effect,  of 
the  placement  of  the  downstream  and  outer  boundaries  will  be  considered 
in  Section  IV. 


Along  the  airfoil  surface,  the  no-slip  adiabatic  conditions 


u  =  V  =  0 


(2.22) 


are  imposed  for  viscous  flows.  For  the  Euler  equations,  the  conditions 


1/  =  0 

Ulh  0  (2-23) 


are  specified,  where  U  ,  ^  are  tJie  contravariant  velocities. 
Equation  (2.16).  The  boundary  conditions  (2.22)  and  (2.23)  are  applied 
in  conjunction  with  a  body-fitted  grid  v^ich  is  nearly  orthogonal  at  the 
airfoil  surface. 

Finally,  along  the  wake-cut  ED,  averaging  is  used  for  all  the  flow 
variables  so  as  to  ensure  continuity. 

Since  only  steady  flows  are  considered  in  this  research,  initial 
conditions  are  not  of  primary  concern  to  provide  a  converged  numerical 
solution.  A  simple  uniform  flow  initial  condition  was  always  employed  in 
the  present  airfoil  calculations,  unless  previous  computed  results  were 
already  available. 


3.  TURBULENCE  MODEL 

Turbulence  is  simulated  by  a  modified  version  of  the  algebraic  eddy 
viscosity  model  of  Baldwin  and  Lomax^.  As  depicted  in  Figure  2,  three 
separate  regions  are  considered  in  the  implementation  of  the  turbulence 
model.  Namely,  the  airfoil  boundary  layers,  the  near-wake  in  the 
vicinity  of  the  airfoil  trailing  edge  and  the  far-wake. 

In  the  airfoil  boundary  layers,  a  two-layer  formulation  is 
employed.  The  inner  turbulent  eddy  viscosity  is  given  by  the 
Prandtl  -  Van  Driest  expression 

e.  =  p(kYD)^l(jl  (2.24) 


D  =  1  -  exp 


1  - 

.  V 

PwKljH  /26 

•  I 

'  J 

(2.25) 


8u  9v 

8y 


(2.26) 


where  w  is  the  vorticity,  Y  represents  the  distance  normal  to  the  airfoil 
surface,  K  =  0.40  is  von  Karman's  constant  and  the  subscript  w  denotes 
val ues  at  the  wall . 

In  the  outer  region  of  the  boundary  layer,  the  turbulent  eddy 
viscosity  is  defined  as  follows 

^0  ~  ^cp  '"max  ^max  ^kleb 


(2.27; 


where  F 


max  (Y|a)|D),  Y 


is  the  value  of  Y  at  which  F 


max 


occurs, 


' 

F 

1  +  5.5 

C  Y/y 

6 

kleb 

kleb  max 

and 


max 


(2.28) 


k  =  0.0168,  C^p  =  1.6. 


^kleb 


0.3 


The  turbulence  model  switches  from  the  inner  to  the  outer 
formulation  at  the  first  value  of  Y  away  from  the  wall  where 
Transition  is  simulated  by  simply  beginning  the  application 
of  the  above  turbulence  model  at  the  boundary  layer  trip  locations 

4 

specified  in  the  experiments. 

The  far-wake  turbulent  eddy  viscosity  is  computed  as  follows 


^wk 


P  c, 


wk 


Y 

^max  ^  DIF 
^max 


'"kleb 


(2.29) 


where 


=  max  (Y|»|)  . 

‘'dIF  =  ♦  ''"’min 


=  0.058 


(2.30) 


and  Y  is  the  value  of  Y  at  which  F  occurs.  In  the  wake,  Y  is 
max  max  ’ 

measured  from  the  wake  centerline  as  determined  from  the  location  of 

minimum  velocity.  The  intermittency  factor  is  obtained  from 

Equation  (2.26)  with  C.,  .  =  0.53.  The  constants  C  ,  and  C.  ,  .  are 

k  1  eb  wk  k  1  eb 

chosen  in  order  to  match  the  above  formulation  with  the  theoretical 
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results  given  by  Schlichting  for  an  incompressible  turbulent  wake 
(see  Appendix  A  for  details). 


The  turbulent  eddy  viscosity  in  the  near-wake  is  computed  by 

allowing  the  trailing  edge  eddy  viscosity  profile  to  exponentially  reach 
its  far-wake  value.  Referring  to  Figure  2,  the  following  expression  is 
employed 

*a[c„,  (x^.Y)  -  e(x,^.Y)] 

*  '  ■  ?)] 

The  distance  x^  -  x^^  is  typically  chosen  to  be  of  the  order  of 
106,  where  6  denotes  the  average  boundary  layer  thickness  at  the  trailing 
edge.  The  sensitivity  of  the  airfoil  lift  and  drag  coefficients  to  the 
variation  of  x^  -  is  addressed  in  Section  IV.  This  near-wake  ad 
hoc  formulation,  similiar  to  that  of  References  14  and  15,  represents  a 
preliminary  approach  within  the  context  of  a  simple  algebraic  eddy 
viscosity  model. 
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SECTION  III 
NUMERICAL  PROCEDURE 


A'.> 

•V^, 


The  numerical  precedure  consists  of  two  basic  steps.  First,  a 
body-fitted  finite-difference  grid,  must  be  constructed  about  the  airfoil 
in  order  to  simplify  the  implementation  of  the  boundary  conditions  and  to 
provide  sufficient  resolution  of  the  flow  features.  Second,  the 

finite-difference  form  of  the  governing  equations  is  solved  using  a 
suitable  numerical  scheme.  In  this  investigation,  grids  were  generated 

O 

by  the  elliptic  technique  of  Visbal  and  Knight.  Two  different  schemes 

were  utilized  for  the  numerical  solution  of  the  flow  equations.  Namely, 

the  implicit  factored  algorithm  of  Beam  and  Warming^^,  and  the  explicit 

g 

unsplit  MacCormack's  algorithm.  The  grid  generation  method,  the 

numerical  schemes  and  the  criterion  for  convergence  of  the  solution  to 

steady  state  are  presented  below. 

1 .  GENERATION  OF  COMPUTATIONAL  GRIDS 

Nearly  orthogonal  body-fitted  grids  were  generated  about  the 

airfoil  using  the  method  developed  by  Visbal^®  and  Visbal  and 

8  17 

Knight.  ’  This  technique,  which  is  based  on  elliptic  partial 
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differential  equations  ,  employs  the  two-step  procedure  described 

below. 

a .  I ntermedi  a  te  T rans  forma ti on 

As  a  first  step,  an  ortiiogonal  grid  is  generated  with  a 
user-prescribed  distribution  of  tiie  mesh  points  along  the  airfoil  and 


s\ 

N* 


>.v. 

V  ' 

*  •, 
V. 

m  *m  ^ 


V  -• 

S’.'*' 

-v* 


wake-cut  (Figure  1).  The  intermediate 
satisfies  the  Poisson  equations^® 

'  5.x  *  5yy  ■  ♦(5-x)  =  j,V 

C  X 


trans forma tion 


[C(x,y),  X(x,y)] 

(3.1) 


V^Y  =  Y  +  Y  “  0 
^  ^  •^xx  -'■yy 

where  , 

hj  =  (x|  +y|)‘ 

\  '  '■'x  "  "x>‘‘ 


(3.2) 


(3.3) 


and  the  forcing  function  (p  (?,x)  results  from  the  general  expression  for 
the  Laplacian  in  orthogonal  curvilinear  coordinates.  In  reference  to 
Figure  1,  the  boundary  conditions  for  ?  and  x  are 


c  = 

■  =  0 
an 

on  Ti 

E  = 

?  0 
^max  ,  an 

on  Fa 

(3.4) 

C  = 

^^(t)  ,  X  =  0 

on  Ts 

ar. 

3n 

■  0  »  X  ■  Xfj,gj( 

on  r4 

where  t 

is  the  arc  length 

along  and  n  denotes  the  normal 

to  the 

corresponding  boundary.  The  boundary  condition  on  simply  indicates 
that  the  grid  points  are  distributed  along  the  airfoil  and  wake-cut  in  a 
desired  monotonic  fashion. 


The  transformation  equations  and  boundary  conditions 
(Equations  (3.1)  -  (3.4))  are  expressed  in  terms  of  (E,x)  derivatives,  as 


shown  in  Reference  17.  The  resulting  equations  are  solved  for 

[x(  5,x)  .y  ( ^.x)]  ill  a  uniform  rectangular  grid  in  the  (5,  x)  plane  using 

point  successive  overrelaxation  (SOR)  and  starting  from  an  arbitrary 

(non-orthogonal )  initial  mesh.  The  forcing  function  <(>  and  the  grid 

points  on  Ij ,  ^  are  dynamically  adjusted  during  the  solution 
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process. 

For  the  present  airfoil  grids,  the  above  intermediate 
transformation  technique  is  applied  only  in  the  vicinity  of  the  airfoil 
(up  to  2.-3.  chords  away).  Since  t^e  grid  is  orthogonal,  it  can  be 
easily  extended  to  the  outer  boundary  using  a  straight-forward 
algebraic  procedure  which  results  in  an  improved  efficiency.  In 
addition,  a  relatively  course  mesh  in  the  x  direction  can  be  utilized, 
b.  Final  Transformation 

In  the  second  step,  a  nearly-orthogonal  mesh  is  constructed 
with  a  user-specified  distribution  of  grid  points  along  T-j  and  T2> 
The  distribution  along  and  is  obtained  from  the  intermediate 
mesh.  Introducing  the  simple  transformation  x  =  X(h),  Equations  (3.1) 
and  (3.2)  become^ ^ 

=  P  (C.n)  =  ^C.x(n)j 


where 


V  =  Q(C.n)  =  (n^  +  n^S 


(3.6) 


IMPLICIT  NAVIER -STOKES  CODE 


The  implicit  code  solves  the  strong  conservative  formulation  of  the 
Navier-Stokes  equations  [Equation  (2.10)]  using  the  approximate- 
factorization  algorithm  of  Beam  and  Warming. This  scheme  in  “delta" 
form  and  with  first-order  Euler  time-differencing  can  be  written  as 
follows  ,  ^ 

fr  r 3a"  fr  .  A.  fae"  3"n]  1  -n 

[  1  1-^  ' 

-At  1^1^  (E,  -  V,  -  V^)"  +  (E^  -  w,  -  Wa)"]  (3.11) 


+  1  n 

I  =  a  + 


(3.12) 


where  n  denotes  the  temporal  index  (i.e.  |  (nAt)),  and 


Jacobian  matrices 


A  =  ^ 


B  = 


(3.13) 


3V 

M  = 

are  given  in  Appendix  B. 


3W 


Application  of  second  order  central  differencing  for  the  space 
derivatives  yields 


I  +  At 


[Vi.j  -  I  ["5  <  =  '  -  ''>>>, J 

■*  >j„  (E,  -  -  «^V,  -  1 
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(3.14) 


(3.15) 


where 


c-  =  (i  -  Da?;  .  1  1  ’  1  IL 


(3.16) 


Oj  =  (j  -  DAn 


1  <  j  <  JL 


Vi.j  ■  •  «n  ^.0  ' 

( 

=  <^n.J  -  w/i.j  =  -  fi.j.ij/Zin 


(3.17) 


are  finite-difference  operators.  The  transformation  derivatives  (x^,  x^, 

y,- .  y  ).  required  in  Equations  (3.14)  and  (3.15),  are  computed  from  the 
5  n 

body-fitted  grid  using  second-order  central  differences  at  interior 
points  and  one-sided  approximations  at  the  boundaries.  The  scheme  is 
implemented  in  a  standard  AOI  fashion  by  solving  Equation  (3.14)  for 
each  n-line  (2^j_<JL-l),  followed  by  the  solution  of  Equation  (3.15)  for 
each  ?-line  (2  £  i  £  IL-1).  This  results  in  the  solution  of  a 
block-tridiagonal  linear  system  along  every  coordinate  line. 

In  order  to  maintain  numerical  stability,  artificial  dissipative 
terms  must  be  added  to  the  basic  Beam-Warming  algorithm. Following 
Reference  19,  explicit  fourth-order  damping 


-0)  AtjT'  .  (6^  +  6D  (U  • 
e  i,j  E  n'  ^i,j 


(3.18) 
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is  appended  to  the  right-hand-side  of  Equation  (3.14).  The  implicit 
second-order  damping  terms 


“(i)  •  At  J 


2 


-(jj^.AtJ 


2 

n 


(3,19) 


are  inserted  within  the  respective  implicit  operators.  The  damping 

coeffi cient^JJ  is  of  order  one  andu  .  >  2^  . 

€  1  —  e 

To  accelerate  convergence  to  steady-state  the  local  time  step 

it.  .  is  employed^® 

* 

(3.20) 


At.  .  =  max 
'  »  J 


At. 


1  + 


'  '^i.j 


where  a  is  constant,  typically  from  1.0  to  2.0.  For  viscous  flows,  At^ 
corresponds  to  a  Courant  number  (CFL)  of  order  50-500  at  the  location  of 
minimum  grid  spacing  on  the  airfoil  (see  Equations  (3.22)-(3.25)  for  the 
definition  of  the  CFL  number). 

Finally,  the  boundary  conditions  are  implemented  in  the  explicit 
fashion  described  in  Reference  12.  A  similar  version  of  the  present 

implicit  code  had  been  previously  validated  for  a  variety  of  flow 

4.-  21 
con  fi  gura  tions. 

3.  EXPLICIT  NAVIER-STOKES  CODE 

The  explicit  Navier-Stokes  code  solves  the  chain-rule  conservative 
form  of  the  governing  equations  (2.1)  utilizing  the  explicit  unsplit 

9 

predictor-corrector  algorithm  of  MacCormack.  Details  of  this 


well-known  algorithm  can  be  found  in  Reference  22  and  therefore  are  not 
included  here.  As  part  of  this  scheme,  a  fourth-order  pressure  damping 
term  is  incorporated  in  order  to  control  numerical  oscillations  in 
regions  of  large  flow  gradients.  This  damping  term,  used  in  both  the 
predictor  and  corrector  steps,  has  the  form 


D  = 

=  -BAtAC' 


(lUl  +  al^l)  T.  I  llpj  33_ 


( 1 1/ 1  +  a  I  vn  1 )  ^1  2lp  1  ia 
3n2  '  3n 


(3.21) 


where  a=fr RT  is  the  speed  of  sound,  U,  U  are  the  contravariant  velocities 
(Equation  (2.16))  and  6  is  the  specified  damping  coefficient  typically 
ranging  from  1.0  to  3.0. 

In  order  to  improve  the  efficiency  for  steady  flow  computations,  a 
local  time  step  is  incorporated  in  the  basic  MacCormack's  algorithm. 
This  is  accomplished  by  specifying  a  constant  Courant  number  (CFL) 
throughout  the  computational  domain.  Namely, 


(x^  +  ytr  AC, 


The  prescribed  Courant  nuirber  is  typically  assigned  a  value  of  0.9. 

Previous  versions  of  this  explicit  Navier-Stokes  code  have  been 

2  23 

successfully  applied  to  a  variety  of  fluid  dynamic  problems.  ’ 
Unlike  the  implicit  code  which  was  written  for  a  scalar  computer,  the 
explicit  solver  is  fully  vectorized  and  exploits  the  vector-processing 
capabilities  of  the  CRAY  computer. 

4.  STEADY-STATE  CONVERGENCE  CRITERION. 

Since  steady  flowfields  are  obtained  by  time-integration  of  the 

Navier-Stokes  equations  from  a  given  initial  condition,  a  suitable 

convergence  criterion  must  be  specified.  Convergence  was  assessed  by 
carefully  monitoring  the  airfoil  lift  (C|^)  and  drag  (C^) 

coefficients.  For  all  cases  computed,  the  monitored  coefficients 
approached  a  steady  value  in  a  damped  oscillatory  fashion.  Here,  the 

flow  was  assumed  converged  when  (1)  the  amplitude  in  the  oscillations  of 
Cj^  was  less  than  0.05  -  0.1%  and  (2)  when  the  variations  in  Cp  were 
within  1-2  drag  counts  (1  drag  count  =  0.0001).  For  the  implicit 
algorithm,  the  corresponding  root -mean- square  values  of  the  residual  were 
typically  10“^  for  all  equations.  Although,  a  more  stringent 
convergence  criterion  could  be  stipulated,  the  present  one  was  found  to 
be  suitable  for  engineering  calculations. 


5. 


CALCULATION  OF  AIRFOIL  FORCE  COEFFICIENTS 

The  airfoil  lift  and  drag  coefficients,  denoted  by  C^  and  C^ 
respectively,  were  computed  by  integration  of  the  pressure  and  shear 
stress  along  the  airfoil  surface.  The  stress  vector  at  the  wall,  fj^, 
is  given  in  cartesian  tensor ial  notation  as  follows 
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k 


n  a, 

1  kl 


(3.26) 


where 


n^  =  Viyivnl 


is  the  normal  to  the  surface  and  is  the  stress  tensor 

0  T 

X  xy 

^xy  *^y 

where  a,=  .peTy^  ''xy  "  yy 

in  Equations  (2.4),  Performing  the  tensorial  inner  product  (3.26)  and 
assuming  u  =  v  =  0  and  e  =  0  at  the  surface,  one  obtains 
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(3.28) 


where  f  ,  f  represent  the  stress  components  in  the  x  and  y 
X  y 

directions  respectively. 


The  total  force  on  the  airfoil  per  unit  span  (F  ,  F  )  is 

X  y 

obtained  by  straight-forward  integration 


using  trapezoidal  rule.  Finally,  the  lift  and  drag  coefficients  are 
given  by 


Cl  =  (-FjjSina  +  F  cosa)/%p^U^ 


(3.29) 


Cp  =  (F^cosa  +  FySina)/isP^U2 


(3.30) 


whereon  is  the  airfoil  angle  of  attack. 


SECTION  IV 


RESULTS  AND  DISCUSSION 


1 .  FLOW  CONFIGURATION  AND  COMPUTATIONAL  DETAILS 

Computations  were  performed  for  an  11%  thick  supercritical  airfoil 

4 

designated  as  DSMA  523.  This  modified  Whitcomb  airfoil  exhibits 
significant  rear-loading  and  a  strong  viscous-inviscid  interaction  as 
compared  to  more  conventional  sections  having  little  or  no  aft-camber. 
This  particular  airfoil  was  selected  as  one  of  the  test  cases  for  the 

5 

1980  Stanford  Conference  on  Complex  Turbulent  Flows.  The  experimental 

4 

database  contains  surface  pressure  distributions  and  velocity  and 
density  profiles  along  the  upper  surface,  lower  surface  and  wake  for  a 
range  of  Mach  nunber,  Reynolds  number  and  angle  of  attack.  Two  specific 


cases 

were  considered  in 

this  study 

.  Namel  y , 

a  subcritical 

case 

(Ko  = 

0.6 ,  a 

=  2.6%  Re^  = 

4  X  10®) 

and  a 

supercritical 

case 

(H»  = 

0.8,  a 

=  1.8',  Re^  = 

2  X  10®). 

The 

corresponding 

flow 

parameters  are 

given  in  detail 

in  Table  1 

• 

In  the 

calculations. 

three  different  computational  grids 

were 

employed.  These  grids,  referred  to  subsequently  as  coarse  (92  x  31), 
medium  (140  x  45)  and  fine  (204  x  55),  were  generated  by  the  procedure 
described  in  Section  III.  As  shown  in  Figure  3  for  the  fine  grid,  the 
mesh  is  nearly  orthogonal  and  displays  substantial  clustering  in  regions 
of  high  flow  gradients.  The  details  of  all  grids  are  summarized  for 
convenience  in  Table  2. 


Numerical  solutions  of  the  flow  equations  were  obtained  utilizing 
either  the  explicit  unsplit  MacCormack 's  scheme  or  the  implicit 

Beam-Warming  algorithm.  Both  schemes  incorporated  a  local  A  t  (Equations 

(3.20)  and  3.22))  in  order  to  accelerate  convergence  to  steady  state. 
The  effect  on  convergence  when  using  a  local  At  was  investigated  for  both 
algorithms  and  is  discussed  next. 

Computations  for  the  airfoil  subcritical  case  on  the  coarse  grid 

using  Beam-Warming  algorithm  indicated  a  gain  in  efficiency  of  about  six 
when  a  local  At  was  employed  (with  At^  =  0.01  ,  o  =  0.5  in  Equation 

(3.20) ).  The  convergence  history  for  the  airfoil  lift  and  drag 

coefficients  is  shown  in  Figure  4  and  points  out  that  an  unphysical 

oscillatory  solution  can  occur  for  higher  values  of  a  and  ^  t^. 

However,  when  a  steady  state  was  reached,  the  results  were  as  expected 

independent  of  a  and  At  . 

0 

MacCormack 's  code  with  a  local  At  displayed  a  speed-up  factor  of 
approximately  seven  for  the  computation  of  a  flat  plate  supersonic 
(Moo  =  3.0)  boundary  layer  with  mass  injection.  Figure  5  shows  the 

convergence  of  the  skin  friction  coefficient  at  a  given  location.  When 
applied  to  the  airfoil,  the  use  of  a  local  At  resulted  in  small  amplitude 
oscillations  about  the  steady  solution  (see  Figure  6).  These 

oscillations  disappeared  when  a  constant  At  was  imposed.  Despite  this 
behaviour,  using  a  local  At  during  the  early  stages  of  the  calculation 
substantially  improved  the  efficiency  of  the  explicit  algorithm. 


TABLE  1 .  FLOW  PARAMETERS 


x/c  tr ans  i  ti  on  Experi  men  tal 


Case 

R^c 

a 

upper 

lower 

®L 

<^0 

Subcri  ti  cal 

0.6 

4x10® 

2.6* 

0.05 

0.18 

0.58 

O.OIl 

Supercritical 

0.8 

1.8* 

0.18 

0.65 

TABLE  2.  GRID  DETAILS 


AS  /c 

AS  /c 

Tl 

No.  of  points 

£ 

Grid 

IL 

X 

JL 

min. 

max. 

min. 

max. 

on  airfoil 

=  0.6 

1 

92 

X 

31 

0.010 

0.060 

0.00025 

1.03 

58 

18.5 

2 

140 

X 

45 

0.005 

0.050 

0.00010 

1.38 

86 

11.5 

3 

204 

X 

55 

0.003 

0.025 

0.00005 

1.25 

150 

6.9 

AS 


2.1 


Legend 

IL,  JL: 


No.  of  grid  points  In  Eand  ndirections  respectively, 
streamwise  grid  spacing  on  airfoil  (Equation  (3.24)) 
grid  spacing  normal  to  airfoil  surface  (Equation  (3.25)) 
normal  spacing  at  airfoil  surface  in  terms  of  law  of 
the  wall  coordinates  (upper  surface,  x/c  =  0.5) 


NOTE:  For  the  above  grids,  the  farfield  boundary  was  located  10  chords 


The  airfoil  calculations  were  performed  on  a  CYBER  175  and  on  the 
NASA  Ames  CRAY  XMP  computer.  The  average  rates  of  data  processing, 
defined  as  CPU  time  per  grid  point  per  iteration,  are  given  in  Table  3 
for  both  algorithms.  The  efficiency  increase  on  the  CRAY  XMP  relative  to 
the  CYBER  175  was  a  factor  of  55  for  the  vectorized  explicit  code. 
However,  the  corresponding  value  for  the  scalar  implicit  algorithm  was 
only  a  factor  of  5.1.  For  a  typical  airfoil  steady  solution  on  a  204  x 
55  grid,  the  implicit  code  was  about  3.8  times  faster  than  the  explicit 
code  on  the  scalar  computer.  However,  the  situation  was  reversed  on  the 
vector  computer,  with  the  explicit  code  being  2.9  times  faster.  It 
should  be  noted,  however,  that  the  efficiency  of  the  implicit  algorithm 
could  be  increased  substantially  by  vectorization  and  by  the  techniques 
of  References  24  and  25. 

A  total  of  23  different  calculations  were  performed  for  the 
selected  airfoil  configuration.  A  detailed  comparison  of  numerical 
results  and  experiments  is  presented  below  for  the  subcritical  and 
supercritical  case  separately. 

TABLE  3. 

COMPARISON  OF  COMPUTER  TIMES  FOR  EXPLICIT 
AND 

IMPLICIT  CODES 


DPRa 

Approx.  No.  of 

Overall  relative  speed^ 

77^ 

A1  gori  thm 

CYBER  175 

CRAY  XMP 

Iterations  req'd 
for  convergence^ 

CYBER  175 

CRAY  XMP 

« 

«  V'  '.*■  * 

BEAM-WARMING 

4.4x10-3  8 

.6x10-^ 

900 

3.8 

0.35 

MACCORMACK 

0.6x10-3  1 

.1x10-5 

25,000 

1 .0 

1.0 

2.  RESULTS  FOR  SUBCRITICAL  CASE 

4 

Since  in  the  experiments  the  ratio  of  wind  tunnel  height  to 
model  chord  was  only  four,  interference  effects  are  expected  to  be 
significant  and  must  be  accounted  for  {assuming  the  data  is  correctable) 
by  adjustments  in  the  freestream  Mach  nuitber  and  angle  of  attack. 
Although  the  main  objective  of  the  present  study  is  the  comparison  of  the 
two  numerical  schemes,  an  attempt  was  made  to  determine  the  effective 
angle  of  attack  for  the  subcritical  case  in  order  to  provide  a  more 
meaningful  comparison  with  the  experiments.  For  this  purpose,  several 
computations  were  performed  with  the  implicit  code  on  the  coarse  grid  at 
various  angles  of  attack.  The  Cp-distribution  and  lift  coefficients, 
shown  in  Figure  7,  indicated  an  effective  angle  of  attack (  a)  of 
approximately  1.7  degrees  in  order  to  match  the  experimental  Cj^.  All 
subsequent  calculations  were  therefore  performed  with  this  value  of  a. 

The  sensitivity  of  the  numerical  solution  to  the  relaxation  length 
xo-Xte  in  the  near  wake  turbulence  model  (Equation  (2.31))  was 
investigated  next.  The  flow  was  computed  with  the  implicit  algorithm 
using  two  values  of  xo-x^^  differing  by  an  order  of  magnitude  (namely, 
0.2c  and  2.0c  where  c  is  the  airfoil  chord  length).  As  Figure  8  shows, 
the  solution  was  insensitive  to  variations  in  the  near-wake  relaxation 
distance.  In  the  remaining  calculations,  the  value  xo-x^p  =  0.2c  was 
therefore  employed.  This  value  is  approximately  86,  where  6  denotes  the 
average  boundary  layer  thickness  at  the  airfoil  trailing  edge. 

According  to  a  formal  truncation  error  analysis,  the  numerical 
smoothing  terms  are  of  high  order  and  therefore  their  effects  on  the 
accuracy  of  the  numerical  solution  are  seldom  considered.  In  the  present 


study,  the  effects  of  artificial  viscosity  on  the  computed  airfoil 
viscous  flowfield  were  investigated  on  a  limited  basis.  Calculations 
were  carried  out  on  the  coarse  mesh  using  both  codes  and  different 
damping  coefficients.  The  results  for  Cp,  C|^  and  are  given  in 
Figure  9.  The  variations  of  Cp  with  damping  were  slightly  larger  for 
Beam-Warming  algorithm  (Figure  9a)  as  compared  to  MacCormack's  scheme 

(Figure  9b).  In  both  cases,  the  lift  coefficient  remained  essentially 
unchanged  while  C^  varied  by  as  much  as  30  drag  counts.  Since  the 

effects  of  damping  are  expected  to  diminish  as  the  mesh  resolution 
increases,  computations  were  performed  on  the  fine  grid  using 

MacCormack's  algorithm  with  two  values  for  the  damping  coefficient 
differing  by  a  factor  of  three.  The  corresponding  variations  in  Cj^  and 
Cjj  with  the  damping  coefficient  were  0.3%  and  eight  counts  respectively. 

In  order  to  investigate  the  effects  of  numerical  resolution,  the 
flowfield  was  computed  on  all  three  grid  systems  (Table  2)  using  both 
Navier-Stokes  codes.  The  results  for  the  implicit  algorithm  are 
displayed  in  Figure  10.  The  canputed  lift  coefficient  increased 
monotonically  with  increasing  resolution  unlike  the  drag  coefficient 
which  did  not  show  any  specific  trend.  The  variations  in  Cj^  and  Cp 
through  the  second  mesh  refinement  were  only  2.2%  and  5  drag  counts 

respectively  as  compared  to  12%  and  19  counts  for  the  first  grid 
refinement.  Although  some  measure  of  numerical  convergence  is  observed, 
the  Cp-distributions  (Figure  10)  indicate  that  further  resolution  on 
the  upper  surface  leading  edge  area  is  still  required.  Figure  11  shows 
the  Mach  number  contours  near  the  relatively  blunt  leading  edge  for  the 
coarse  and  fine  grids.  Besides  being  much  more  smooth,  the  fine  grid 


solution  makes  apparent  the  existence  of  an  embedded  supersonic  region 
(terminated  by  a  shock)  at  the  leading  edge.  The  corresponding  mesh 
refinement  results  for  the  explicit  scheme,  shown  in  Figures  12  and  13, 
displayed  a  behaviour  similar  to  that  of  the  implicit  algorithm.  For  the 
explicit  code,  the  variations  in  and  Cp  after  the  second  grid 
refinement  were  1.0%  and  15  drag  counts  respectively. 

A  comparison  between  MacCormack  and  Beam-Warming  results  as  well  as 
between  the  computed  results  and  the  experiments  is  presented  below. 
Figures  14,  15  and  16  contain  the  results  for  Cp,  C|^  and  Cp  on  the 
coarse,  medium  and  fine  grids  respectively.  The  agreement  between  the 
two  calculated  Cp-distributions  improves  as  the  mesh  is  refined. 
However,  differences  still  persist  on  the  airfoil  upper  surface 
immediately  downstream  of  the  eiibedded  supersonic  region.  The  computed 
lift  and  drag  coefficients  obtained  with  the  two  algorithms  on  the  fine 
grid  (Figure  16)  differ  by  1.8%  and  5  drag  counts  respectively.  These 
discrepancies  are  acceptable  for  engineering  applications.  The  better 
agreement  (in  terms  of  and  Cp)  between  the  two  numerical  schemes 
on  the  coarse  mesh  (Figure  14)  can  only  be  regarded  as  fortuitous. 

In  order  to  illustrate  the  importance  of  viscous  effects  for 
supercritical  airfoils,  the  inviscid  flowfield  was  computed  with  the 
implicit  code  on  the  coarse  (Figure  14)  and  medium  (Figure  15)  grids. 
Even  for  a  small  angle  of  attack,  subcritical  flow  without  significant 
boundary  layer  separation,  the  decrease  in  due  to  viscous  effects  is 
substantial  (approximately  29%  for  the  medium-grid  solution.  Figure  15) 

Referring  to  Figure  15,  the  discrepancies  in  Cp  and  Cj^  between 
the  fine-grid  computational  results  and  the  experiments  are  in  part  due 


to  uncertainties  in  the  effective  freestream  Mach  number  and  angle  of 
attack.  The  better  prediction  of  the  experimental  lift  coefficient  on 
the  coarse  mesh  (Figure  14)  simply  points  out  that  comparison  with 
experimental  data  can  be  misleading  unless  proper  numerical  resolution 
criteria  are  first  established. 

A  detailed  comparison  of  computed  and  experimental  velocity 
profiles  is  given  in  Figures  17  -  19  for  the  airfoil  upper  surface,  lower 
surface  and  wake.  As  expected,  the  nimerical  profiles  obtained  with  the 
explicit  and  implicit  codes  are  in  very  good  agreement  with  each  other  at 
all  locations.  On  the  suction  surface  (Figure  17),  the  predicted 
profiles  show  higher  values  for  the  freestream  velocity.  This  is 
consistent  with  the  higher  computed  (Figure  16)  and  can  be 

attributed  mainly  to  uncertainties  in  the  effective  angle  of  attack.  The 
discrepancies  between  the  measured  and  computed  velocity  profiles  at  the 
airfoil  cove  (Figure  18)  and  in  the  near-wake  (Figure  19)  are  most  likely 
due  to  deficiencies  in  the  algebraic  turbulent  eddy  viscosity  model 
employed. 

A  qualitative  comparison  of  computed  density  contours  and  the 
experimental  interferogram  is  presented  in  Figure  20.  All  typical 

features  of  the  subcritical  flow  about  the  airfoil  are  discernible  and  in 
good  qualitative  agreement  with  the  experimental,  flowfield.  A  close-up 
comparison  near  the  trailing  edge  is  shown  in  Figure  21. 

It  is  interesting  to  compare  the  uncertainties  in  the  numerical 

simulation  of  different  algorithms  with  the  experimental  uncertainties 
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associated  with  different  wind  tunnel  occupancy  periods.  Figure  22a 
shows  two  measured  C„-distributions  for  the  same  nominal  subcritical 


conditions.  The  specific  comparison  indicates  that  the  experimental 
uncertainty  is  of  the  same  order  of  magnitude  as  the  discrepancies 

between  the  two  numerical  algorithms  (Figure  16). 

3.  RESULTS  FOR  SUPERCRITICAL  CASE 

Calculations  for  the  supercritical  flow  case  (Table  1)  using  both 
numerical  algorithms  were  limited  to  the  fine  grid  (Table  2).  The 
coarser  grid  systems  are  not  expected  to  provide  sufficient  spatial 
resolution. 

In  order  to  achieve  a  more  meaningful  comparison  between 

computations  and  experiments,  an  attempt  was  made  to  determine  the 
effective  frees tream  Mach  number  and  angle  of  attack.  As  shown  in 
Figure  23,  computations  at  the  nominal  conditions  {Moo=  0.8,  a  =1.8") 
predicted  a  shock  location  much  further  aft  as  compared  to  the 

experimental  data.  The  freestream  Mach  number  M<»=  .765  and  angle  of 

attack  a  =  0.7*  were  found  to  approximate  the  experimental  shock 

position.  These  values  were  therefore  employed  in  all  subsequent 

calculations.  As  suggested  in  Reference  5,  the  experimental  Cp  and 
were  corrected  to  account  for  the  shift  in  Mach  nunber  by  keeping  a 
constant  static -to -total  pressure  ratio.  The  present  computations 

indicated  that  this  supercritical  airfoil  is  extremely  sensitive  to  small 
variations  in  the  freestream  Mach  number. 

Since  the  farfield  boundary  conditions  employed  in  this  study  are 
only  approximate  (see  Section  II),  the  effect  of  the  outer  boundary 

placement  on  the  numerical  solution  must  be  investigated.  For  this 
purpose,  the  supercritical  flowfield  was  computed  with  the  explicit 
algorithm  for  three  different  locations  of  the  farfield  boundary  (namely, 
10,  25  and  50  chords  away).  The  results  for  C„,  C,  and  are 
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shown  in  Figure  24  and  indicate  that  as  the  outer  boundary  is  placed 

further  avay  from  the  airfoil  the  shock  moves  downstream  with  a 

corresponding  increase  in  lift.  The  difference  in  C|^  between  the  25 

chord  and  50  chord  solution  is  only  1.8%.  Although  additional  outward 
boundary  placements  have  not  been  considered,  it  may  be  concluded  that  a 
computational  domain  extending  at  least  50  chords  away  from  the  airfoil 
should  be  employed  for  strong  viscous-inviscid  interacting  calculations. 
The  farfield  boundary  effects  might  be  less  significant  for  more 

conventional  airfoils  or  for  improved  formulations  of  the  boundary 
condi  ti  ons . 

A  comparison  of  the  computed  results  and  experiments  is  presented 
below.  Figure  25  displays  reasonable  agreement  between  the  two  numerical 
Cp-distributions  except  in  the  vicinity  of  the  shock  where  additional 
spatial  resolution  in  the  streamwise  direction  is  still  needed.  Both 
computed  solutions  failed  to  predict  the  measured  static  pressure  on  the 
upper  surface  downstream  of  the  shock  and  on  the  lower  surface  cove.  It 
appears  that  this  discrepancies  point  out  deficiencies  in  the  algebraic 
turbulence  model  downstream  of  the  shock-boundary-layer  interaction  and 
in  the  cantered  region.  A  comparison  of  tJie  velocity  profiles  on  the 
upper  surface,  lower  surface  and  wake  is  given  in  Figures  26,  27  and  28 
respectively.  The  computed  velocities  are,  as  expected,  in  close 
agreement  with  each  other.  Discrepancies  between  the  experimental  and 
computed  velocity  profiles  in  the  cove  and  near-wake  are  apparent. 
Finally,  a  qualitative  comparison  of  the  experimental  and  numerical 
flowfields  is  presented  in  Figure  29  in  terms  of  computed  density 
contours  and  the  experimental  interferogram. 


It  is  interesting  to  compare  the  uncertainties  in  the  numerical 
results  of  different  algorithms  with  the  repeatiabil  i  ty  in  the 
experimental  measurements.  Figure  22b,  reprinted  from  Reference  4,  shows 
the  pressure  distributions  obtained  for  the  same  nominal  conditions 

during  two  different  wind  tunnel  runs.  For  this  case,  the  experimental 
uncertainty  is  found  to  be  of  the  same  order  of  magnitude  as  that 
encountered  in  the  numerical  simulation  (Figure  25). 

The  thin-layer  approximation  (Equation  (2.20))  has  been  commonly 
employed  in  the  numerical  solution  of  high  Reynolds  nunber  airfoil 

or 

flows,  and  justification  for  its  use  can  be  found  in  Reference  7. 
For  unsteady  transonic  airfoils,  significant  discrepancies  between  the 
thin-layer  and  the  full  Navier-Stokes  results  were  observed  by  Chyu  and 

p  o  pft 

Kuwahara.  On  the  other  hand,  Degani  and  Steger  found  good 

agreement  between  the  two  formulations  when  applied  to  a  supersonic 

compression  ramp.  In  the  present  study,  the  effect  of  the  thin-layer 

approxi nation  was  investigated  for  the  airfoil  supercritical  flow.  The 

results,  shown  in  Figure  31,  indicated  very  good  agreement  between  the 

tliin-layer  and  the  full  Navier-Stokes  solution.  Therefore,  the  use  of 
the  thin-layer  equations  is  fully  justified  for  this  particular  type  of 
fl  ows . 

To  illustrate  the  extent  of  the  viscous-inviscid  interaction 
effects  for  aft-canbered  airfoils,  the  flowfield  was  also  computed  using 
the  Euler  equations.  As  Figure  31  shows,  the  presence  of  the  boundary 
layer  modifies  the  flowfield  dramatically.  Viscous  effects  result  in  a 
much  weaker  shock  (displaced  upstream)  and  in  a  substantial  decrease  in 

the  "plateau"  loading.  This  is  a  direct  consequence  of  boundary  layer 


value 


the  continuity  equation  are  given  in  Table  4  (the  results  for  the 
remaining  equations  were  found  to  be  similar).  The  root -mean- square 
value  of  the  damping  relative  magnitude  was  always  less  than  one. 
However,  the  maximum  value  was  an  order  of  magnitude  greater.  A  contour 
plot  of  the  damping-to-truncation  error  ratio  (Equation  (4.2))  for  the 
continuity  equation  is  shown  in  Figure  31  for  the  supercritical  case. 
Only  contour  values  greater  than  one  were  drawn  for  clarity.  As 
expected,  the  region  of  large  error  is  located  in  the  vicinity  of  the 
shock,  where  the  accuracy  of  the  solution  is  already  degraded. 

A  second  assessment  of  the  numerical  error  due  to  damping  was 
performed  by  comparing  the  artificial  viscosity,  introduced  by  the 
smoothing  terms,  with  the  real  physical  viscosity  y  +  e.  This  was  done 
for  the  explicit  algorithm  due  to  the  particular  form  of  the  damping 
terms  (Equation  (3.21)).  Comparison  of  the  damping  term  with  the 
normal  diffusive  terms  in  the  streamwise  momentum  equation  (Equations 
(2.1)  -  2.4))  yields  the  following  approximate  expression  for  the 

artificial-to-real  viscosity  ratio. 
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Contours  of  the  above  ratio  are  shown  in  Figure  32  for  both  the 
subcritical  and  supercritical  case.  It  can  be  seen  than  within  the 
airfoil  viscous  region  the  artificial  viscosity  is  one  order  of  magnitude 
smaller  than  the  true  viscosity  y  +  e  .  A  similar  analysis  can  be 
performed  for  the  C -direction  damping  term  D^.  This  term  was  only  found 
to  be  significant  for  the  supercritical  case  in  the  vicinity  of  the  shock 


where  becomes  comparable  in  magnitude  to  the  normal  diffusive  term. 


The  above  results  can  only  be  regarded  as  a  preliminary  attempt  to 
assess,  a  posteriori,  the  accuracy  of  the  computed  flowfields,  and 
further  grid  refinement  still  represents  the  only  conclusive  approach  for 
determining  appropriate  numerical  resolution  criteria. 


TABLE  4.  DAMP I NG-TO-TR UNCATION  ERROR  RATIO 


FOR  THE  CONTINUITY  EQUATION 


SECTION  V 


CONCLUSIONS  AND  RECOMMENDATIONS 


In  the  present  study,  a  critical  examination  of  the  numerical 
solution  of  high  Reynolds  number  transonic  airfoils  was  performed. 
Solutions  of  the  flowfield  about  a  supercritical  airfoil  were  generated 
by  solving  the  two-dimensional  Navier-Stokes  equations  with  an  algebraic 
turbulent  eddy  viscosity  model.  The  governing  equations  were  solved  on 
curvilinear  body-fitted  grids  using  two  different  numerical  algorithms. 
Namely,  the  explicit  unsplit  MacCormack's  scheme  and  the  implicit 
Beam-Warming  algorithm.  Both  subcritical  and  supercritical  flows  were 
considered,  with  a  total  of  23  cases  computed. 

The  numerical  uncertainties  associated  with  different  schemes,  grid 
resolution,  artificial  viscosity  and  farfield  boundary  placement  were 
carefully  examined.  Numerical  simulations  produced  by  the  two  numerical 
schemes  with  different  conservative  formulations  exhibited  nearly 
identical  velocity  profiles  and  density  contours.  The  maximum 
discrepancy  between  these  solutions  was  generally  less  than  the 
repeatable  data  scattering  band  of  experimental  measurements.  However, 
the  uncertainty  in  the  more  stringent  requirement  of  the  computed  lift 
and  drag  coefficient  is  still  of  the  order  of  2..0  and  20  drag  counts 
respectively  for  the  finest  grid  (204  x  55)  utilized.  Further  numerical 
resolution,  particularly  in  shock  regions,  is  still  required.  This  could 
perhaps  be  achieved,  in  a  more  efficient  manner,  by  the  use  of  grid 
adaptation  techniques. 


/II 


A  comparison  of  the  supercritical  flowfields  obtained  with  the 
Euler,  thin-layer  and  Navier -Stokes  equations  was  performed.  The 
thin-layer  approximation  gave  essentially  the  same  results  as  the  full 
Navier-Stokes  equations  and  therefore  its  use,  for  this  type  of  flows,  is 
clearly  justified.  However,  numerical  solutions  of  the  Euler  equations 
failed  to  provide  a  reasonable  approximation  of  the  flowfield  due  to  the 
dramatic  viscous/inviscid  interaction  effects  for  transonic  aft-cambered 
airfoils . 

Comparison  of  the  computed  results  and  the  experimental  data  showed 
a  reasonable  prediction  of  all  of  the  essential  features  of  the  flow. 
However,  detailed  comparison  of  velocity  profiles  on  the  airfoil  upper 
surface  lower  surface  and  wake  pointed  out  deficiencies  of  the  turbulence 
model  in  the  cove  region,  downstream  of  the  shock/boundary  layer 
interaction  and  in  the  near-wake. 

The  comparison  of  computed  and  experimental  results  was  clouded  by 
uncertainties  in  both  the  measurements  and  in  the  numerical  procedure. 
In  principle,  the  numerical  uncertainties  could  be  reduced  to  a  desired 
acceptable  level.  However,  the  uncertainties  in  the  measurements  were 
not  adequately  documented  in  the  experimental  study.  For  the 
investigated  supercritical  airfoil,  the  problem  is  significant  due  to  the 
extreme  sensitivity  of  the  flow  to  variations  in  the  freestream  Mach 
number  and  angle  of  attack.  Therefore,  before  a  conclusive  evaluation  of 
turbulence  modeling  can  be  conducted,  an  experimental  study  aimed  at  the 
determination  of  interference  effects  should  be  performed.  An 
alternative  and  less  expensive  approach  of  evaluating  interference 
effects  would  be  the  numerical  simulation  of  the  entire  three-dimensional 


wind-tunnel  flowfield. 


APPENDIX  A 
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CONSTANTS  FOR  FAR-WAKE  TURBULENCE  MODEL 


The  similarity  velocity  profile  in  a  two-dimensional  incompressible 


turbulent  wake  can  be  expressed  as  follows 


13 
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where 


n  .  1/2  1/2  y 


Un-c  =  U  -  U  .  =  6.3535  ( - 7- 

DiF  max  mm  px 


1/2 


is  the  turbulent  eddy  viscosity  and  (=Uoo)  and  are  the 

maximum  (freestream)  and  minimum  velocities,  respectively. 


Assuming 


9U/3Y 


in  Equation  (2.30),  and  using  Equation 


(A.l)  for  the  velocity,  it  can  be  easily  shown  that 


F  =  2UQ.p  ^  exp  (-n^) 
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^wk  0.0579  ^0 

Therefore,  in  order  to  have  the  constand  must  be 

set  equal  to  0.0579. 

For  the  wake,  the  constant  appearing  in  Equation  (2.28)  is 

determined  as  follows 
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where  b  is  the  wake  half-width 
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and  Y  is  given  in  Equation  (A.3).  The  resulting  value  for  C.  . 
max 

is  0.527. 


APPENDI X  B 
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JACOBIAN  MATRICES  FOR  THE  NAVIER-STOKES  EQUATIONS 


The  Jacobian  matrices  defined  in  Equation  (3.13)  are  given  below. 
These  matrices  are  obtained  by  straight-forward  differentiation,  after 
rewriting  the  vectors  E-j,  E2,  W2  (Equations  (2.11)  -  (2.15)) 

in  terms  of  the  conserved  variables  (  p/J,  pU/J,  pv/J,  pe/J).  The 
Jacobian  matrix  A  is 


where 
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and  U  is  given  in  Equation  (2.16).  The  Jacobian  matrix  B  can  be  obtained 
from  Equation  (B.l  )  by  applying  the  substitutions 
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The  Jacobian  matrix  M  is 
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where  the  coefficients  b^.  (i  =1,  .  .  .,4)  are  given  in  Equation 
(2.17).  The  Jacobian  matrix  N  can  be  obtained  from  Equation  (B.2)  by 
replacing  the  coefficients  b^.  with  the  corresponding  coefficients  d^. 
(Equation  (2.19)). 
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FIGURE  2  Regions  for  Turbulence  Model 
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Convergence  of  Skin  Friction  Coefficient  for  Explicit  Code 
(Supersonic  Flat  Plate  Boundary  Layer  with  Blowing) 
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FIGURE  6 


Convergence  of  Ci  and  Co  for  Explicit  Code 
(Airfoil  Subcritical  Case,  Fine  Grid) 


Resolution  on  Computed  Subcritical  Flow,  Implicit  Code 
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FIGURE  15  Computed  and  Experimental  Cp,  Ci  and  Cn  for  Subcritical  Case 
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FIGURE  25  Computed 
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